#include "../src/analyzer_common.h"

int main(int argc, char **argv)
{
    int op = 0;
    char mtx_path[100];
    if (argc >= 3)
    {
        strcpy(mtx_path, argv[1]);
        op = atoi(argv[2]);
    }
    else if (argc >= 2)
    {
        strcpy(mtx_path, argv[1]);
    }
    else
    {
        printf("%s /path/to/matrix [op_type]\n", argv[0]);
        exit(1);
    }

    sparse_matrix_t *A = init_smatrix();
    sparse_matrix_t *U = init_smatrix();
    sparse_matrix_t *L = init_smatrix();
    read_matrix(mtx_path, A);

    struct timeval time_start;
    double time;
    sparse_value_t *rhs = NULL;
    sparse_matrix_t *A_prune = NULL;
    sparse_index_t n;
    sparse_value_t *sol, *rhs_chk, *rhs_tmp;
    sparse_index_t *diagptr;
    sparse_matrix_t *A_bsr;
    int block_order;

    switch (op)
    {
    case 0: // matrix analyze
        matrix_analyze(A);
        break;
    case 1:                       // export upper
        get_triangle(1, 0, A, U); // upper
        sprintf(buf, "%s_U.lid", mtx_path);
        export_lid(buf, U);
        break;
    case 2:                       // export lower
        get_triangle(0, 0, A, L); // lower
        sprintf(buf, "%s_L.lid", mtx_path);
        export_lid(buf, L);
        break;
    case 3: // convert to lid
        sprintf(buf, "%s.lid", mtx_path);
        export_lid(buf, A);
        break;
    case 4: // convert to mtx
        sprintf(buf, "%s.mtx", mtx_path);
        export_mtx(buf, A);
        break;
    case 5: // triangle matrix level-set detect
        level_set_detect(A);
        break;
    case 6: // submtx [0:argv[3])
        sprintf(buf, "%s_0_%d.mtx", mtx_path, atoi(argv[3]));
        analyzer_submtx(A, 0, atoi(argv[3]), buf);
        break;
    case 7:
        analyzer_droprate(A, atoi(argv[3]));
        break;
    case 8:
        sprintf(buf, "%s_block_%d.mtx", mtx_path, atoi(argv[3]));
        analyzer_submtx_blockfirst(A, atoi(argv[3]), buf);
        break;
    case 9:
        if (argc >= 4)
        {
            sparse_index_t n;
            read_vector(argv[3], &n, &rhs);

            n /= 3;
            sparse_value_t *rhs_sub = malloc(sizeof(sparse_value_t) * n);
            for (sparse_index_t row = 0; row < n; row++)
            {
                rhs_sub[row] = rhs[row * 3];
            }
            free(rhs);
            rhs = rhs_sub;
        }
        else
        {
            rhs = malloc(sizeof(sparse_value_t) * A->m);
            for (sparse_index_t row = 0; row < A->m; row++)
            {
                rhs[row] = 0;
                for (sparse_pointer_t idx = A->rowpointer[row]; idx < A->rowpointer[row + 1]; idx++)
                {
                    rhs[row] += A->value[idx];
                }
            }
        }
        analyzer_diagsolver(A, rhs);
        break;
    case 10:
        A_prune = malloc(sizeof(sparse_matrix_t));
        memset(A_prune, 0, sizeof(sparse_matrix_t));
        analyzer_prune(A, A_prune);
        sprintf(buf, "%s_prune.lid", mtx_path);
        printf("prune done. Expoting to %s\n", buf);
        export_lid(buf, A_prune);
        break;
    case 11:
        read_vector(argv[3], &n, &sol);
        read_vector(argv[4], &n, &rhs_tmp);
        n /= 3; // block_order == 3
        rhs_chk = malloc(sizeof(sparse_value_t) * n);
        for (sparse_index_t row = 0; row < n; row++)
        {
            rhs_chk[row] = rhs_tmp[row * 3];
        }
        free(rhs_tmp);
        analyzer_diagsolver_check(A, sol, rhs_chk);
        free(rhs_chk);
        free(sol);
        break;
    case 12:
        diagptr = malloc(sizeof(sparse_index_t) * A->m);
        analyzer_timer_start(&time_start);
        #pragma omp parallel for
        for (sparse_index_t i = 0; i < A->m; ++i)
        {
            for (sparse_pointer_t j = A->rowpointer[i]; j < A->rowpointer[i + 1]; ++j)
            {
                if (A->columnindex[j] == i)
                {
                    diagptr[i] = j;
                    break;
                }
            }
        }
        time = analyzer_timer_stop(&time_start);
        printf("ilu0 setup time : %.2lf ms\n", time * 1e3);

        analyzer_timer_start(&time_start);
        ilu0_ref(A->rowpointer, A->columnindex, A->value, A->m, A->nnz, diagptr);
        time = analyzer_timer_stop(&time_start);
        printf("ilu0 factorise time : %.2lf ms\n", time * 1e3);

        free(diagptr);

        // sprintf(buf, "%s_ilu0.lid", mtx_path);
        // export_lid(buf, A);
        break;
    case 13:
        analyzer_blas_check();
        break;
    case 14: // bsr ilu0
        block_order = atoi(argv[3]);
        A_bsr = malloc(sizeof(sparse_matrix_t));
        analyzer_timer_start(&time_start);
        analyzer_csr_to_bsr(A, A_bsr, block_order);
        time = analyzer_timer_stop(&time_start);
        printf("bilu0 csr->bsr time : %.2lf ms\n", time * 1e3);

        free(A->rowpointer);
        free(A->columnindex);
        free(A->value);
        free(A);
        A = NULL;

        analyzer_timer_start(&time_start);
        #pragma omp parallel for
        for(sparse_pointer_t inidx = 0; inidx < A_bsr->nnz * block_order * block_order; inidx++){
            A_bsr->value[inidx] += 0.0;
        }
        time = analyzer_timer_stop(&time_start);
        printf("bilu0 memtouch time : %.2lf ms\n", time * 1e3);

        diagptr = malloc(sizeof(sparse_index_t) * A_bsr->m);
        analyzer_timer_start(&time_start);
        #pragma omp parallel for
        for (sparse_index_t i = 0; i < A_bsr->m; ++i)
        {
            for (sparse_pointer_t j = A_bsr->rowpointer[i]; j < A_bsr->rowpointer[i + 1]; ++j)
            {
                if (A_bsr->columnindex[j] == i)
                {
                    diagptr[i] = j;
                    break;
                }
            }
        }
        time = analyzer_timer_stop(&time_start);
        printf("bilu0 setup time : %.2lf ms\n", time * 1e3);

        analyzer_timer_start(&time_start);
        bilu0_ref(A_bsr->rowpointer, A_bsr->columnindex, A_bsr->value, A_bsr->m, A_bsr->nnz, diagptr, atoi(argv[3]));
        time = analyzer_timer_stop(&time_start);
        printf("bilu0 factorise time : %.2lf ms\n", time * 1e3);

        free(diagptr);
        break;
    case 15:
        block_order = atoi(argv[3]);
        analyzer_structed_sparse(A, block_order);
        break;
    }
}
